# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import functools
from hysop.tools.htypes import check_instance, first_not_None
from hysop.tools.decorators import debug
from hysop.tools.numpywrappers import npw
from hysop.backend.host.host_operator import HostOperator, OpenClMappable
from hysop.core.graph.graph import op_apply
from hysop.operator.base.solenoidal_projection import SolenoidalProjectionOperatorBase
from hysop import __DEFAULT_NUMBA_TARGET__
[docs]
def build_projection_filter(FIN, FOUT, K, KK, target=__DEFAULT_NUMBA_TARGET__):
assert len(FIN) == len(FOUT) == 3
assert len(K) == len(KK) == 9
args = FIN + K + KK + FOUT
try:
import numba as nb
from hysop.tools.numba_utils import make_numba_signature, prange
signature, _ = make_numba_signature(*args)
layout = "(m,n,p),(m,n,p),(m,n,p), "
layout += "(m),(n),(p), (m),(n),(p), (m),(n),(p), "
layout += "(m),(n),(p), (m),(n),(p), (m),(n),(p) "
layout += "-> (m,n,p),(m,n,p),(m,n,p)"
@nb.guvectorize([signature], layout, target=target, nopython=True, cache=True)
def filter_projection_3d(
Fin0,
Fin1,
Fin2,
K00,
K01,
K02,
K10,
K11,
K12,
K20,
K21,
K22,
KK00,
KK01,
KK02,
KK10,
KK11,
KK12,
KK20,
KK21,
KK22,
Fout0,
Fout1,
Fout2,
):
for i in prange(0, Fin0.shape[0]):
for j in prange(0, Fin0.shape[1]):
for k in range(0, Fin0.shape[2]):
F0 = Fin0[i, j, k]
F1 = Fin1[i, j, k]
F2 = Fin2[i, j, k]
G0 = K00[i] * F0
G1 = K11[j] * F1
G2 = K22[k] * F2
L0 = KK00[i] * F0
L1 = KK11[j] * F1
L2 = KK22[k] * F2
if (i != 0) or (j != 0) or (k != 0):
C0 = (L0 + K10[i] * G1 + K20[i] * G2) / (
KK00[i] + KK01[j] + KK02[k]
)
C1 = (K01[j] * G0 + L1 + K21[j] * G2) / (
KK10[i] + KK11[j] + KK12[k]
)
C2 = (K02[k] * G0 + K12[k] * G1 + L2) / (
KK20[i] + KK21[j] + KK22[k]
)
else:
C0 = 0
C1 = 0
C2 = 0
Fout0[i, j, k] = Fin0[i, j, k] - C0
Fout1[i, j, k] = Fin1[i, j, k] - C1
Fout2[i, j, k] = Fin2[i, j, k] - C2
return functools.partial(filter_projection_3d, *args)
except ModuleNotFoundError:
def filter_projection_3d(
Fin0,
Fin1,
Fin2,
K00,
K01,
K02,
K10,
K11,
K12,
K20,
K21,
K22,
KK00,
KK01,
KK02,
KK10,
KK11,
KK12,
KK20,
KK21,
KK22,
Fout0,
Fout1,
Fout2,
):
for i in range(0, Fin0.shape[0]):
for j in range(0, Fin0.shape[1]):
for k in range(0, Fin0.shape[2]):
F0 = Fin0[i, j, k]
F1 = Fin1[i, j, k]
F2 = Fin2[i, j, k]
G0 = K00[i] * F0
G1 = K11[j] * F1
G2 = K22[k] * F2
L0 = KK00[i] * F0
L1 = KK11[j] * F1
L2 = KK22[k] * F2
if (i != 0) or (j != 0) or (k != 0):
C0 = (L0 + K10[i] * G1 + K20[i] * G2) / (
KK00[i] + KK01[j] + KK02[k]
)
C1 = (K01[j] * G0 + L1 + K21[j] * G2) / (
KK10[i] + KK11[j] + KK12[k]
)
C2 = (K02[k] * G0 + K12[k] * G1 + L2) / (
KK20[i] + KK21[j] + KK22[k]
)
else:
C0 = 0
C1 = 0
C2 = 0
Fout0[i, j, k] = Fin0[i, j, k] - C0
Fout1[i, j, k] = Fin1[i, j, k] - C1
Fout2[i, j, k] = Fin2[i, j, k] - C2
return functools.partial(filter_projection_3d, *args)
[docs]
def make_div_filter(target, *args):
try:
import numba as nb
from hysop.tools.numba_utils import make_numba_signature, prange
from hysop import __DEFAULT_NUMBA_TARGET__
signature, _ = make_numba_signature(*args)
layout = "(m,n,p), (m,n,p), (m,n,p), (m),(n),(p) -> (m,n,p)"
@nb.guvectorize([signature], layout, target=target, nopython=True, cache=True)
def compute_div_3d(Fin0, Fin1, Fin2, K0, K1, K2, Fout):
for i in prange(0, Fin0.shape[0]):
for j in prange(0, Fin0.shape[1]):
for k in range(0, Fin0.shape[2]):
Fout[i, j, k] = (
K0[i] * Fin0[i, j, k]
+ K1[j] * Fin1[i, j, k]
+ K2[k] * Fin2[i, j, k]
)
except ModuleNotFoundError:
def compute_div_3d(Fin0, Fin1, Fin2, K0, K1, K2, Fout):
for i in range(0, Fin0.shape[0]):
for j in range(0, Fin0.shape[1]):
for k in range(0, Fin0.shape[2]):
Fout[i, j, k] = (
K0[i] * Fin0[i, j, k]
+ K1[j] * Fin1[i, j, k]
+ K2[k] * Fin2[i, j, k]
)
return functools.partial(compute_div_3d, *args)
[docs]
class PythonSolenoidalProjection(
SolenoidalProjectionOperatorBase, OpenClMappable, HostOperator
):
"""
Solves solenoidal projection (project a 3d field F such that div(F)=0),
using spectral methods.
"""
[docs]
def build_divergence_filters(self, target=__DEFAULT_NUMBA_TARGET__):
if self.compute_divFin:
FIN, K, DIV_IN = self.FIN, self.K, self.DIV_IN
K = (K[0], K[4], K[8])
args = FIN + K + DIV_IN
self.pre_filter_div = make_div_filter(target, *args)
else:
self.pre_filter_div = None
if self.compute_divFout:
FOUT, K, DIV_OUT = self.FOUT, self.K, self.DIV_OUT
K = (K[0], K[4], K[8])
args = FOUT + K + DIV_OUT
self.post_filter_div = make_div_filter(target, *args)
else:
self.post_filter_div = None
[docs]
@debug
def setup(self, work):
super().setup(work=work)
FIN, FOUT = self.FIN, self.FOUT
K, KK = self.K, self.KK
self.filter_projection = build_projection_filter(FIN, FOUT, K, KK)
self.build_divergence_filters()
@op_apply
def apply(self, simulation=None, **kwds):
super().apply(**kwds)
with self.map_objects_to_host():
for Ft in self.forward_transforms:
Ft()
if self.compute_divFin:
self.pre_filter_div()
self.backward_divFin_transform()
self.filter_projection()
if self.compute_divFout:
self.post_filter_div()
self.backward_divFout_transform()
for Bt in self.backward_transforms:
Bt()
if self.compute_divFin:
self.ddivFin.exchange_ghosts()
if self.compute_divFout:
self.ddivFout.exchange_ghosts()
for Bt in self.backward_transforms:
Bt.dfield.exchange_ghosts()